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ABSTRACT 

We revisit the computation of the extragalactic gamraa-ray signal from cosmological 
dark matter annihilations. The prediction of this signal is notoriously model depen- 
dent, due to different descriptions of the dumpiness of the dark matter distribution 
at small scales, responsible for an enhancement with respect to the smoothly dis- 
tributed case. We show how a direct computation of this "flux multiplier" in terms 
of the nonlinear power spectrum offers a conceptually simpler approach and may ease 
some problems, such as the extrapolation issue. In fact very simple analytical recipes 
to construct the power spectrum yield results similar to the popular Halo Model 
expectations, with a straightforward alternative estimate of errors. For this specific 
application, one also obviates to the need of identifying (often literature-dependent) 
concepts entering the Halo Model, to compare different simulations. 

Key words: Cosmology: theory-dark matter-gamma-rays: diffuse background 



1 INTRODUCTION 

To unveil the particle physics nature of the Dark Matter 
(DM) inferred by astrophysical and cosmological observa- 
tions, a plethora of strategies are currently pursued, includ- 
ing production in high energy colliders and direct detection 
of recoils in underground detectors (see the monograph in 
|Bertone||2010| for a review). However, indirect detection of 
Galactic and extragalactic by-products of DM annihilations 
(or decays) is the sole way to access information on the 
DM remotely, i.e. in the environments where evidence for 
its presence has been collected. In particular, with the DM 
in the form of Weakly Interacting Particles (WIMPs) signals 
in gamma-rays are expected from the inner halo or center of 
the Milky Way or from Milky Way halo substructures (ei- 
ther known via their baryonic counterparts, dwarf galaxies, 
or as "dark satellites"). All these signals, while promising, 
depend on the "local" (in cosmological terms) DM environ- 
ment, related for example to the assembly history of our 
Galaxy. On the other hand, at high Galactic latitudes the 
diffuse signal from the Milky Way halo is expected to be 
roughly comparable with the diffuse extragalactic one. This 
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Extragalactic DM annihilation Flux (EDMF) depends in a 
statistical sense on how DM is distributed at different scales, 
and how this distribution evolves with redshift. Although 
typically more challenging to detect than its Galactic coun- 
terpart, the EDMF retains important cosmological informa- 
tion otherwise impossible to access. Since the underlying 
particle physics parameters are the same, one might hope to 
corroborate an eventual signal of DM annihilation by look- 
ing at this high-latitude flux, or to isolate the cosmological 
information by comparison with the Galactic signal. 

The purpose of this article is to revisit the computation 
of the "flux multiplier" , which is the main poorly known in- 
gredient relating the EDMF to the parameters determining 
the Galactic flux. Although this issue has been subject to a 



lot of attention in different contexts (see e.g. Bergstrom et al 
2001[|Ullio et al.|2002[|Silk fc Taylor|2003[ ) here we advocate 
a more straightforward computation via direct integration 
of the nonlinear power spectrum (PS). Although formally 
equivalent to the Halo Model (HM) framework in conflgura- 
tion space, especially when discussing individual component 
uncertainties and cosmology-dependency, this approach has 
the great advantage of involving mostly model-independent 
quantities (i.e. whose definitions do not rely on the specific 
parameterization of the variables and quantities entering the 
HM) which can be extracted directly from simulations. In 
addition, some uncertainties on clustering properties are eas- 
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ier to assess in the this context. Even when a HM approach 
is used to guide the extrapolation to the very small scales, 
the PS approach has the benefit to isolate the "effective 
combination of parameters" responsible for the signal. 

This article is structured as follows. In Sec. [2] we intro- 
duce the central notion of flux multiplier ({z) and recap the 
main formulae used for its computation within the HM. In 
Sec.[3]we show how the "systematic" uncertainty in the sig- 
nal can be more effectively reformulated in terms of the PS 
extrapolation. Finally, in Sec. |4] we discuss our results and 
comment on possible future directions. 



of each halo (halo density profile and concentration). This 
allows one to write (^{z) as in |Ullio et al.| p002| 

f dM^M^^{F) , (5) 

with dn/dM the comoving density of halos per unit mass, 
A„ the mean halo over density, and F being the function 

1 , , fn" dxx^K^(x) , , 

FiM,z)^ct{M,z) /l ' (6) 



2 STANDARD FORMALISM 

Formally, for a constant annihilation cross section (av) the 
flux (number of photons per energy interval, unit area, time 
and solid angle) from DM particles with mass can be 
written (see e.g. |Ando fc Komatsu|2006| ) 

-r(..E) p2(^^^{l-^ dN{E',z'^ 



4>(E, n) = 



(av) c 
8 TV mj. 



I 



dz- 



H{z) (14-2)3 



dE' 



(1) 



with H[z) the Hubble expansion function, p referring to the 
DM density at redshift z in the angular direction Cl, e~^'^'^^ 
accounting for absorption onto the extragalactic background 
light, and dN/dE being the emitted spectrum per annihila- 
tion, with E' denoting the energy corresponding to present- 
day value E at redshift z. In general this spectrum includes 
not only the prompt photons, from lines or tt" decays, but 
also secondary emission from energy losses of other particles 
(such as inverse-compton from onto the CMB). Equa- 
tion ([T]) is usually rewritten isolating the directional depen- 
dence contained in the density contrast field 5 as: 

p{z, d) = [l + 5{z, d)] p{z) = [1 + S{z, Cl)]noM Pc{l + zf . (2) 

Due to the deeply nonlinear regime of matter perturbations 
nowadays, a very good approximation is to neglect the "av- 
erage" matter contribution in the equation above (the "1" 
term in square brackets) and concentrate on the one due to 
the dumpiness (the S term). In addition, if, as here, the an- 
gular dependence is not of interest and all one cares about 
is the average fiux over the angular direction, one has 

where we defined the flux multiplier 

az) ^ {s'\z,ci)) , (4) 

namely the variance of density fiuctuations over the sky at a 
given epoch. The largest uncertainty in the EDMF compu- 
tation stems from (^{z), on which we concentrate henceforth. 
In order to compute the above quantity, one customarily 



resorts to the Halo Model (HM) framework (see Cooray & 
Sheth 2002 for a review) . The HM assumes that all the mass 
in the Universe is contained in virialized objects (halos) fully 
characterized by their mass. As a consequence, the statistical 
properties of the mass density field are determined by the 
spatial distribution of matter within an individual halo and 
by the spatial distribution of halos, assumed not to overlap 
one with respect to the other. For this calculation the crucial 
quantities are given by the number density of halos of a given 
mass [halo mass function) and by the density distribution 



which depends on the DM halo density p = K{M, z)Hi{x), 
where K includes the cosmology dependence in terms of 
the variables {M, z} and k{x) is the assumed universal 
shape function determined by numerical simulations. Also, 
one assumes that the halo density is non-vanishing only 
within a radius Rv (virial radius), which is conventionally 
parametrized via the concentration parameter c„ = Rv/rs\ 
the dimensionless variable x is just r/r,,. The halo properties 
need not to be universal: in that case F has to be intended 
as an average over the probability distribution of the rele- 
vant parameters (most notably c^). Finally, is commonly 
normalized by imposing that all mass resides in some halo. 

In general, a more faithful description of simulations 
requires accounting for halo "sub-structures" (sub-halos) 
which in turn have their own mass-function, concentration 
and shape/profile properties. Also, one can distinguish dif- 
ferent "halo populations" , according to the degree of dynam- 
ical interaction they undergo with other structures: the ones 
the HM is typically compared to are the so-called "distinct" 
or "isolated" halos. All in all, in order to perform an estimate 
of C,{z), the usual practice is to fit from simulations the fol- 
lowing quantities: ^^{M,z), A^{z), Cy{M,z) (as well as its 
distribution around the mean), k{x) (as well as its distribu- 
tion around the mean and possible 2— evolution), and similar 
quantities for each different category of objects: sub-halos, 
non-distinct halos, etc. Additionally, some approximations 
(e.g. spherical shape of the halos) are implicitly assumed. 



3 DIRECT COMPUTATION OF ( FROM THE 
POWER SPECTRUM 

While the HM has proven extremely useful for the semi- 
analytical modeling of nonlinear clustering, the EDMF does 
not really depend directly on the many different parame- 
ters and functions of the HM, but only on the PS at very 
small scales. In fact, the flux multiplier of Eq. (j4]| can also 
be written as the r — >■ limit of the two-point correlation 
function (2PCF) = (S {x + r) S (x)). Perhaps this is best seen 
if one thinks of fluctuations of the DM field in momentum 
space and rewrites ("(z) in terms of the nonlinear matter 
power spectrum Pnl — defined as a Fourier transform of the 
2PCF— as 



lim 

7 i-0 



d k sin kr 
k kr 



-Ak,z). 



(7) 
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Figure 1. Dimensionless power spectrum A(fc) [data points) from the Millennium Simulation jSpringel et al.|2005| l at z = {left panel, 
blue) and at z = 2 {right panel, red). The dotted curve denotes the corresponding linear power spectrum, the continuous curve the 
halof it ( [Smith et al.|[2003| prediction while the three different EIM predictions of GBSC are given by the long-dashed curves for the 
CI model, the dot-dashed for C2 and the short-dashed for C3. The gray, straight line represent the the shot-noise limit expected if the 
simulation particles were a Poisson sampling from a smooth underlying density field. 



Here Ajvi:,(fc) is the dimensionless nonlinear PS defined as[^ 
AjvL(fe) = k'^ Pn L{k) / {2n^), while the effective ultraviolet 
(UV) cutoff kmax correspouds to the acoustic oscillation or 
free-streaming damping scale kcut, dependent on the DM 
candidate's kinetic decoupling temperature from the plasma 
in the early Universe. Typical expectations for WIMPs are 
in the range kcut G [10^, 10*] h Mpc"'^ at the end of the linear 
regime of structure formation (see |Green et al.|2005| [Bring- 



2009 



and refs. therein), which are usually mapped 
into corresponding "smallest halo" masses of the order of 
Mmi^ ~ 47r/3(7r/fc^ax)^ G [10-l^ 10"^] MqF] In the follow- 
ing, we shall consider kmax ~ 10^/iMpc~ as a represen- 
tative value. Notice that the low-fc region contributes quite 
insignificantly to the above integral, since the fluctuations 
of the field are still in the perturbative regime. The size of 
the observable Universe provides the physical infrared cut- 
off, fcnim ~ 0(10"^) /iMpc-\ 

Unfortunately, there is no rigorous theoretical way 
to compute AjvL(fc, z). Linear theory can be extended in 
perturbation theory only to the mildly nonlinear regime 
(fe < 0.3 /i Mpc"'^ at 2; = 0) where the matter density field 
can be treated as a fluid in the single-stream approxima- 
tion (see Bernardeau et aL]|2002 for a review). At smaller 



Note that our convention for Piv 



■M 



is linked to the 



P^j^(k) defined e.g. in (Bernardeau et al. 2002 i via PNL{k) 



4 TrPj^^ (fc)/(27r)^; the A]vi(fc)'s are instead the same and corre- 
spond to standard definitions. 

In the HM, UV cutoffs enter both via the r -> limit of the 
halo profile (HP) and via the minimum halo mass. When ex- 
trapolating HPs found in current numerical simulations one finds 
that the inner halo (e.g. r <Si scale radius rs in the popular NFW 
HP ( [Navarro et al.|1996[ l) contributes only modestly to the global 
annihilation signal; thus Mmin de facto regulates the extrapola- 
tion uncertainty. Also note that while the mapping between a 
cutoff in fc-space and M„iin ha^s some ambiguities that are intrin- 
sic to the HM, the associated error is subleading compared to the 
extrapolation error. 



scales, the occurrence of shell-crossing and the eventual viri- 
alization of over-dense regions prevents the problem to be 
treated rigorously from first principles, and simulations are 
necessary for quantitative inferences. Yet, theoretical con- 
siderations based on the theory of spherical collapse and 
stable clustering, that is the idea that virialized regions de- 
couple from Hubble expansion keeping a constant physical 
size, allowed [Hamilton et al.| ( |l991| ) to derive a mapping be- 
tween the linear and the nonlinear matter correlation func- 
tion, eventually extended to the PS in [Peacock fc Dod^ 
(|1994[|1996| ). Similarly, the popular halof it formula ( [Smith 
et al.|2003 1 provides a cosmology-dependent mapping (fitted 



to N-body simulations) from the linear to nonlinear PS, sig- 
nificantly improving the accuracy of the earlier approaches. 

Of course, one can also extrapolate the nonlinear PS 
predicted within the HM approach. In this framework, the 
PS can be decomposed into the sum of a term coming from 
the correlation of two DM particles within the same halo 
{1-halo term) and a term describing the correlation between 
particles belonging to distinct halos {2-halo term), related to 
the two-point function of the halo distribution. The highly 
nonlinear regime is clearly described by the 1-halo contri- 
bution which dominates the enhancement factor (. Notice 
that the 1-halo contribution depends on the same ingredi- 
ents required by the usual calculation of the annihilation flux 
described in the previous section, just expressed in fc-space. 
If the HM could reproduce completely the PS extracted from 
simulations, going through the HM or a direct computation 
via integration of the PS would produce identical results. 
However, it is known that several aspects are missed in the 
HM, especially in describing deeply nonlinear scales. For in- 
stance [Giocoli et al.[ ( [2010[ hereafter GBSC) extends the 
basic picture to include the effect of a subhalo population 
and the scatter of the concentration parameter. 

It is thus useful to compute ({z) by performing the in- 
tegral in Eq. ([7|, relying as much as possible directly on 
the results of simulations. A typical result for the nonlin- 
ear matter PS measured in simulations at 2 = and at 
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Table 1. Enhancement factor (divided by 10"^) at four different redshifts from the Millennium Simulation (MS) power spectrum 
and its extrapolations, the halofit formula (HF), and three models of GBSC , denoted CI, C2, C3 as in the original work. 
Results from the extrapolations in Eq.s ijsjl and ^ are given in the form "min"/"max", respectively . Columns 2 to 6 assume an 
integration up to /cmax = 50/iMpc~^; columns (7 to 11) assume /cmax — 10'^ h Mpc ^ while columns 12 and 13 present the result 
of extending the integration up to fcmax = 10^ hMpc~^ . Column 13 is like column 12, but for the Millennium-II Simulations. 



z = 2, is reproduced by the data points in Fig. [TJ adapted 
from Springel et al. (20051. This is compared with the lin- 



ear PS prediction for the same cosmology (dotted curves) 
and the nonlinear predictions of the halofit fitting formula 
(continuous curves). Table |3] shows the enhancement factor 
^ computed using the halofit prediction and assuming a 
sharp cut-ofF in the integral given by kmax = 50, 10* and 
10^ Mq in the columns 3, 8, and 13, respectively and 
for 2; = 0, 1 and 2. For instance, we obtain — 2.7 x 10'^ for 
kmax ~ 50/iMpc~^ at 2 = 0, very close to the simulation 
result^ Lacking a firm theory, we can tentatively assume 
the above extrapolations as "fiducial" prediction for the PS. 
We also show three other HM predictions corresponding to 
different mass-concentration relations and sub-halo models, 
denoted CI, C2 and C3 as in GBSC , with the long-dashed, 
dot-dashed and short-dashed curves, respectively. In addi- 
tion, the gray straight line shows the resolution limit of the 
simulations due to the shot-noise contribution. The values 
of ^ calculated in terms of the simulations results assuming 
kmax = 50/iMpc~^ are given in column 2 of Table [s] for 
z = 0, 1 and 2, showing a good agreement among numerical 
results and fits/models. But, clearly, the behavior of the DM 
field at scales well below the resolution presently achieved 
by N-body simulations is dominating the fiux. In order to 
estimate the error from an extrapolation of the simulations 
results we can make a very simple assumption inspired by 
direct inspection of the data: AjvL(fc) is a non-decreasing 
function of k, but its second derivative is negative at suffi- 
ciently large scales. Hence, we assume that the true Ajv_l(A;) 
is bracketed by the two following cases 

A^r(fe) = AML(k*) , forfc > k* , (8) 
A^r(fc) = AiVL(fc*)+AVt(fc*)(fc-fc*) , for > fc* , (9) 

where k* denotes the scale beyond which the spectrum can- 
not be trusted anymore because of numerical resolution ef- 
fects, and the prime stands for derivative with respect to 
k. For illustration, in columns 7 and 12 of Table [S] we re- 
port in the format "min"/"max" the results of such ex- 
trapolations done beyond k* = {29,16,9.8} /iMpc~^ for 
z — 0,1,2, respectively. The fc*'s have been fixed here as 
the values where the expected shot-noise component reaches 
1% of ANL(k* , z). The extrapolation is then performed cor- 
recting the data for shot-noise. The columns 7 and 12 of 

^ In our calculation we add to the halofit prescription the cor- 
rection, relevant at small scales, provided by J. Peacock here: 
http://www.roe.ac. uk/~jap/haloes/| 



Table [3] show how ^(2) depends in fact on the kmax value. 
Note that, although the functional form of the fitting func- 
tions of Eq. ([8| and Eq. ([9| is heuristical, several different 
approaches to modeling nonlinear gravitational clustering 
suggest a fiattening of A(fc) at large k, i.e. A'(fc) — >■ at 



large k (see e.ff. |Padmanabhan fc Ray 2006 Angulo & White 
2010 1, giving some confidence that Eq. (|8j) and Eq. (|9j) rep- 



resent extreme cases. 

The results in Table [3] illustrate the following points: 

• At least up to fcmax = 10*/i/Mpc (middle set of 
columns), the comparison between the the halofit formula 
and the extended HM is reasonable, with halofit expecta- 
tions falling in the middle of the GBSC predictions at low 
z, while slightly below at higher z\ notice however that HM 
predictions are typically overshooting simulations at high z 
(see e.g. Klypin et al. 20111. Additionally, the simple ex- 
trapolations of Eq.s (|8| and ([9| nicely bracket these model 
expectations and confirms a posteriori that they can be seen 
as reasonable error estimates on the signal. 

• The last columns represent our best guess for a conser- 
vative error range and best estimate of C,(z), assuming the 
representative particle physics input kmax = 10^/iMpc~^. 
Based on the MS, the value ^(0) « lO'^'* appears as a con- 
servative upper limit, and C(Q) ~ lO**'^ a conservative lower- 
limit. On the other hand, a similar preliminary analysis for 
the result of the Millennium-II simulation extracted from 



Fig. 6 in ( Boylan-Kolchin et al. 2009 1 yields one order of 



magnitude reduction in the extrapolation error, as shown in 
column 13, thanks to the higher resolution (we estimated 
k* = {212,120,84} /iMpc"^ for z = 0,1,2). Less conserva- 
tive estimates of k* or a more careful treatment may further 
reduce the error. 

Of course, our results do not include effects unaccounted 
for in the pure cold DM models considered here, nor addi- 
tional particle physics and cosmology assumptions. Without 
the need to discuss any uncertainty on auxiliary variables 
such as concentration, mass function, substructure, etc. we 
reproduced the range of uncertainty computed with tradi- 



tional methods in configuration space (see e.g. Abdo et al. 
|2009[ Fig. 1) with the MS data, while improving by one or- 
der of magnitude over it when using MS-II data. Also the 
decreasing functional shape of ^(2) is in rough agreement 
with existing computations. Obviously, this is a simpler and 
more direct method than what currently employed in the 
literature and it makes obvious that no information on sep- 
arate ingredients of the HM can be inferred by EDMF. 
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4 DISCUSSION AND CONCLUSIONS 

In the context of WIMP models for DM, the extragalactic 
flux from the residual annihilation of these particles repre- 
sents one possible handle for the discovery of the particle 
physics nature of DM, as well as perhaps the only tool to 
extract cosmological informations from very small scales. In 
this article we revisited the issues involved in the compu- 
tation of the main cosmological ingredient for this signal, 
the so-called flux multiplier ({z). Our results suggest that a 
direct computation of ({z) in terms of nonlinear power spec- 
trum presents several advantages with respect to traditional 
HM computations in configuration space. In short, it allows 
one to 

• isolate the essential physical quantity which the signal 
depends on, merely ANL{k, z), rather than requiring to fit 
several theoretical functions, which in the HM parameterize 
halo-related observables which are useful in many instances, 
but irrelevant for the present application; 

• bypass the problems linked to ill-defined concepts or 
different algorithms used by different research groups and 
publications to isolate each auxiliary variable: the PS con- 
stitute in principle a complete description of nonlinear mat- 
ter clustering relevant for ({z) (of course, bearing in mind 
the problem of knowing the UV behaviour, see discussion in 
Sec. [3] and footnote 2); 

• ease the extrapolation problem: one is reduced to find a 
plausible extrapolation for a single function and for "only" 
about 4 decades in k (under the conventional assumptions 
on the UV behaviour of ANL{k,z) referred to above) mak- 
ing the estimate much less prone to cumulative systematic 
effects of a wrong choice of several fitting functions. Even 
when working with the HM, one should gauge the plausi- 
bility of different extrapolations by comparison with the PS 
yielded by rmmerical simulations; to put it simply such in- 
formation IS available and should not be discarded; 

• ease the task of estimating a credible error budget for 
Ciz): simple extrapolations of ANLik,z), suggested by in- 
spection of its behavior in the resolved regime, proved to 
be as spread as (for MS) and significantly better (for MS- 
II) than the existing range of HM estimates. This promising 
results will of course require more extensive analyses and 
inspection of a wider set of simulations (as well as for effects 
like the ones of baryons, neutrinos, etc.), which we leave for 
future work. 

Also, one might envisage to develop further insights 
along related directions. One example would be to assess 
the extent to which one may combine information from 
simulations available for structures at different scales, e.g. 
high-resolution cosmological simulation such as Millennium- 



II and galactic-scales simulations such as Aquarius ( Springel 



et al. 20081 and Via Lactea (Kuhlen et al. 20081. Note 



that the cosmological PS is already commonly obtained 
from simulations by separating "large scales" and "smaller 
scales" ( Jenkins et al.|1998 Springel et al.|2005 l: by choos- 
ing the small cells so that they contain statistically the same 
initial (linear) power, each one contributes about the same 
to the final power, which allows for a relatively sparse sam- 
pling. Similarly, as long as one carries out high-resolution 
simulations of "representative" cells it might be possible to 
extend the range of the nonlinear PS which can be directly 
extracted from simulations. 



Additionally, one might extend the reasoning to com- 
pute the anisotropy in the EDMF: since its 2PCF is related 
to the variance of the DM power spectrum, computing the 
"Ci's" only requires additionally the "squeezed" four point 
correlation function. Although some early attempts to de- 
rive the signal purely on the basis of simulation results ex- 
ist ( Cuoco et al.|[2007 l, it would be interesting to combine 
some recent theoretical insights with modern simulations. 
In fact, compared with the average EDMF, a larger fraction 
of the measurable anisotropy signal should follow from the 
weakly nonlinear scales (see e.ff. |Ando fc Komatsu| ( [2006[ ).) 

Finally, although the HM has proven quite flexible and 
effective in describing most scales of interest for cosmology in 
the last decade, there are some theoretical arguments why 
it might be less and less suitable to describe the large— fc 
and high-z behavior of AjvL(fc,2). As argued for example 
in Elahi et al. (20091, at smaller and smaller scales haloes 



do not fully virialize before being accreted, with the bound- 
ary of (sub) structures becoming increasingly ill defined. Ex- 
trapolations based on mass functions computed at galactic 
scales to the bottom of the CDM hierarchy are thus ques- 
tionable. It would be important to decouple the expected 
growing uncertainties and inadequacies of a specific (albeit 
till now successful) model from the actual dependence of the 
signal of interest here from the underlying cosmology, which 
is more neatly encoded in a direct computation in terms of 
the power spectrum or two-point correlation function. 
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